//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2008 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : Cg4Solver.h
//----------------------------------------------------------------------------
// EW: This solver is superseded by RicSolver, because taking omega=1 reduces Ric to Mic=Cg4.

#ifndef CG4SOLVER_H
#define CG4SOLVER_H

#include "defines.h"
#include "Solver.h"
#include "Array2D.h"

#define CG4SOLVER_BASE Solver

/// conjugate gradient solver
class Cg4Solver : public CG4SOLVER_BASE
{
    typedef CG4SOLVER_BASE base;

public:

    Cg4Solver(Array2D<REAL> *pmatC,
              Array2D<REAL> *pmatS,
              Array2D<REAL> *pmatW,
              int N,
              int M,
              REAL tolerance);

    virtual ~Cg4Solver();

    virtual void solve(Array2D<REAL> *pB,
                       Array2D<REAL> *pX);

private:

    /// @cond nevercall
    Cg4Solver(const Cg4Solver& m);
    Cg4Solver & operator =(const Cg4Solver& m);
    /// @endcond

    void preconditioner();

    Array2D<REAL> m_a1; ///< preconditioned array matW
    Array2D<REAL> m_a2; ///< preconditioned array matS
    Array2D<REAL> m_a3; ///< preconditioned array matC
    Array2D<REAL> m_D; ///< temporary array
    Array2D<REAL> m_Y; ///< temporary array
    Array2D<REAL> m_Z; ///< temporary array
    Array2D<REAL> m_Az; ///< temporary array

};

inline Solver* NewCg4Solver(Array2D<REAL> *pmatC, ///< central matrix
                            Array2D<REAL> *pmatS, ///< south matrix
                            Array2D<REAL> *pmatW, ///< west matrix
                            int N, ///< width
                            int M, ///< height
                            REAL tolerance) ///< solver tolerance
{
    return new Cg4Solver(pmatC, pmatS, pmatW, N, M, tolerance);
}

#endif // CG4SOLVER_H
